This study uses landslide, fire, and precipitation data to study the antecedent precipitation of post-wildfire landslides. The impact of wildfire on rainfall-triggered landslide susceptibility in a variety of global regions is evaluated through a proxy of normalized precipitation.
Use this section to modify file paths if running the notebook on a different computer.
# Landslide data from the NASA Global Landslide Catalog
slide.path <- '../regions_data/glc/GLC20180821.csv'
# Fire data from the MODIS Burned Area dataset
modis.paths <- list.files(path='../regions_data/processed',
pattern='slide_modisburndate_[ns][ew].csv',
full.names = T)
# Precipitation data from CHIRPS
precip.chirps.ww.path <- file.path("../regions_data/processed",
"glc_precip_global.csv")
precip.chirps.ww.monthly.path <- file.path("../regions_data/", "processed",
"glc_precip_global_monthly.csv")
# Additional preciptiation data from Daymet for a limited region
precip.daymet.ca.path <- file.path("../regions_data/processed",
"glc_precip_daymet.csv")
Include landslides with the following parameters: * Rainfall triggered
rain.triggers <- c('rain', 'downpour', 'flooding', 'continuous_rain')
start.date <- lubridate::ymd('2000-01-01')
end.date <- lubridate::ymd('2019-01-01')
min.latitude <- -50
max.latitude <- 50
min.location.accuracy <- '10km'
modis.start.date <- lubridate::ymd('20040101')
modis.end.date <- lubridate::ymd('20190101')
wilcox.alternative <- 'less' # n - burned, m - unburned
wilcox.alpha <- 0.05
# Identify if a location is in the global limits
is.in.global.study <- function (x, y) { between(y, min.latitude, max.latitude) }
# Identify if a location is in the CA/NV pilot study
get.data('state.boundaries', tigris::states)
boundaries <- subset(state.boundaries, STUSPS %in% c('CA', 'NV'))
boundaries.df <- fortify(boundaries)
is.in.ca.study <- function (x, y) {
bbox <- sf::st_bbox(subset(state.boundaries, STUSPS == 'CA'))
between(x, bbox$xmin, bbox$xmax) &
between(y, bbox$ymin, bbox$ymax)
}
Seasons are split in the Northern Hemisphere by:
In the southern hemisphere, the split is shifted instead:
A set of functions converts months and days-of-the-year to seasons.
season.breaks <- c('spring'=59, 'summer'=151, 'fall'=243, 'winter'=334)
yday.to.season <- function(yday, latitude) {
if (latitude > 0) {
case_when(
yday > season.breaks['winter'] ~ 'Winter',
yday > season.breaks['fall'] ~ 'Fall',
yday > season.breaks['summer'] ~ 'Summer',
yday > season.breaks['spring'] ~ 'Spring',
TRUE ~ 'Winter')
} else {
case_when(
yday > season.breaks['winter'] ~ 'Summer',
yday > season.breaks['fall'] ~ 'Spring',
yday > season.breaks['summer'] ~ 'Winter',
yday > season.breaks['spring'] ~ 'Fall',
TRUE ~ 'Summer')
}
}
month.to.season <- function (month, north) {
if (north) {
case_when(
month >= 12 ~ 'Winter',
month >= 9 ~ 'Fall',
month >= 6 ~ 'Summer',
month >= 3 ~ 'Spring',
TRUE ~ 'Winter')
} else {
case_when(
month >= 12 ~ 'Summer',
month >= 9 ~ 'Spring',
month >= 6 ~ 'Winter',
month >= 3 ~ 'Fall',
TRUE ~ 'Summer')
}
}
# Region colors
region.plot.order <- c('All landslides',
'California',
'Intermountain West',
'Pacific Northwest',
'Central America',
'Himalayas',
'Southeast Asia')
scale_region_values <- c('All landslides' = 'grey50',
setNames(RColorBrewer::brewer.pal(6, 'Dark2'),
region.plot.order[2:7]))
scale_region_labels <- c('All landslides' = 'None',
setNames(region.plot.order[2:7],
region.plot.order[2:7]))
scale_color_region <- scale_color_manual('Region',
values = scale_region_values,
labels = scale_region_labels)
scale_fill_region <- scale_fill_manual('Region',
values = scale_region_values,
labels = scale_region_labels)
scale_color_precip <- scale_color_distiller('Precipitation percentile',
palette='BrBG', direction = 1,
guide = guide_colorbar(barwidth = unit(0.5, 'npc')))
scale_fill_precip <- scale_fill_distiller('Precipitation percentile',
palette='BrBG', direction = 1,
guide = guide_colorbar(barwidth = unit(0.5, 'npc')))
# Fire colors
burned.labels <- c('burned'='Burned', 'unburned'='Unburned')
burned.colors <- c('burned' = '#d95f02', unburned = '#7570b3')
scale_color_burned <- scale_color_manual('',
values = burned.colors,
labels = burned.labels)
scale_fill_burned <- scale_fill_manual('',
values = burned.colors,
labels = burned.labels)
# Season
scale_x_season <- scale_x_continuous(
breaks = season.breaks,
labels = c('Spring', 'Summer', 'Fall', 'Winter'))
col.ocean <- 'lightskyblue1'
theme_map <- list(
theme_bw(),
coord_quickmap(),
labs(x='', y='')
)
theme_global_map <- list(theme_map,
borders(fill = 'grey90', size = 0.1),
theme(panel.background = element_rect(fill=col.ocean),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.margin = unit(c(0,0,0,0), "null")),
coord_quickmap(ylim = c(-60, 75), xlim = c(-180, 180))
)
theme_panel <- theme(panel.background = element_rect(fill=col.ocean),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.margin = unit(c(1,0,-1,0), "lines"),
legend.position = 'null')
theme_panel_legend <- theme(
legend.box.margin = margin(-1, 0, 1, 0),
legend.box.just = c(0, 0),
#legend.background = element_rect(color='black'), #debug
legend.position = 'bottom'
)
capitalize <- function(string) {
substr(string, 1, 1) <- toupper(substr(string, 1, 1))
string
}
load.slides <- function (slide.path,
start.date, end.date,
is.in.global.study,
rain.triggers,
min.location.accuracy) {
# Define factor levels
ls.sizes <- c("small", "medium", "large", "very_large", "catastrophic")
ls.accuracy <- c("exact", "1km", "5km", "10km",
"25km", "50km", "100km", "250km")
# Read data file and remove duplicates
slide.raw <- readr::read_csv(slide.path,
guess_max = 100000,
na = c("", "unknown", "--"),
col_types = cols(
landslide_size = readr::col_factor(ls.sizes, ordered=TRUE),
location_accuracy = readr::col_factor(ls.accuracy, ordered=TRUE),
landslide_trigger = readr::col_factor(NULL),
landslide_category = readr::col_factor(NULL),
landslide_setting = readr::col_factor(NULL),
country_code = readr::col_factor(NULL)
)
) %>%
dplyr::mutate(event_date = as_date(ymd_hm(event_date)),
submitted_date = as_date(ymd_hm(submitted_date)),
last_edited_date = as_date(ymd_hm(last_edited_date)))
slide.unique <- slide.raw %>%
dplyr::group_by(event_date, latitude, longitude) %>%
dplyr::mutate(i = row_number()) %>%
dplyr::filter(i == 1) %>%
dplyr::select(-i) %>%
dplyr::ungroup()
# Filter landslides to meet study parameters
slide.unique %>%
dplyr::filter(event_date >= start.date,
event_date <= end.date,
is.in.global.study(longitude, latitude),
landslide_trigger %in% rain.triggers,
location_accuracy <= min.location.accuracy) %>%
droplevels() %>%
dplyr::mutate(event_yday = yday(event_date),
debris_flow = landslide_category %in% c('mudslide', 'debris flow'),
season = factor(map2_chr(event_yday, latitude,
~yday.to.season(.x, .y)),
levels = c('Spring', 'Summer', 'Fall', 'Winter'),
ordered = T))
}
get.data('slide.all.precip', overwrite = F,
load.slides,
slide.path, start.date, end.date,
is.in.global.study, rain.triggers,
min.location.accuracy)
slide.all.precip %>%
dplyr::select(landslide_trigger, location_accuracy, latitude, event_date) %>%
summary()
## landslide_trigger location_accuracy latitude
## rain :1823 exact: 466 Min. :-46.77
## continuous_rain: 588 1km :1707 1st Qu.: 10.49
## downpour :3254 5km :2437 Median : 29.91
## flooding : 48 10km :1103 Mean : 23.52
## 3rd Qu.: 39.11
## Max. : 49.98
## event_date
## Min. :2000-10-13
## 1st Qu.:2010-11-30
## Median :2013-06-24
## Mean :2013-05-17
## 3rd Qu.:2016-03-11
## Max. :2018-08-07
load.precip.all <- function (precip.chirps.ww.path) {
precip.chirps.ww.raw <- read_csv(precip.chirps.ww.path,
col_types = cols(X1 = col_datetime(format = "%Y-%m-%d")))
precip.chirps.ww.raw %>%
dplyr::select(date=X1, OBJECTID, pctl=precip_pctl_7) %>%
dplyr::filter(OBJECTID %in% slide$OBJECTID) %>%
dplyr::mutate(window=7) %>%
dplyr::filter(!is.na(pctl)) # Remove NA values from beginning of data range
}
get.data('precip.chirps.ww.all', precip.chirps.ww.path)
precip.chirps.ww.all
date.filter.precip <- function (raw, slide=slide.all.precip, pre=90, post=30) {
raw %>%
dplyr::right_join(slide %>% dplyr::select(OBJECTID, event_date),
by='OBJECTID') %>%
dplyr::select(OBJECTID, dplyr::starts_with('precip_'), event_date, date) %>%
dplyr::mutate(date = lubridate::as_date(date),
days_before = as.integer(event_date - date)) %>%
dplyr::filter(days_before < pre, days_before > -post)
}
tidy.precip <- function (precip) {
precip %>%
dplyr::select(OBJECTID, starts_with('precip_'), event_date, date) %>%
gather(key='variable', value='amount', starts_with('precip_')) %>%
separate(variable, into=c(NA, 'variable', 'window')) %>%
dplyr::filter(variable=='pctl') %>%
dplyr::select(-variable, pctl=amount) %>%
dplyr::mutate(days_before = as.integer(event_date - date),
yday = yday(event_date - days(days_before)),
window = as.numeric(window)) %>%
dplyr::select(-event_date)
}
load.event.precip <- function (precip.chirps.ww.path) {
precip.chirps.ww.raw <- read_csv(precip.chirps.ww.path,
col_types = readr::cols(X1 = col_datetime(format = "%Y-%m-%d")))
precip.chirps.ww.raw.filt <- precip.chirps.ww.raw %>%
dplyr::rename(date = X1) %>%
dplyr::filter(year(date) > 2003)
precip.chirps.ww.raw.filt2 <- precip.chirps.ww.raw.filt %>%
date.filter.precip
precip.chirps.ww.raw.filt2 %>%
tidy.precip
}
get.data('precip.chirps.ww.event', precip.chirps.ww.path)
precip.chirps.ww.event
load.precip.monthly <- function (monthly.path, slide) {
read_csv(precip.chirps.ww.monthly.path,
col_types = cols(X1 = col_date(format = '%Y-%m-%d'))) %>%
dplyr::transmute(OBJECTID,
year = year(X1), month = month(X1),
precip.mm=precip) %>%
dplyr::inner_join(slide %>% select(OBJECTID, latitude), by='OBJECTID') %>%
dplyr::mutate(north = latitude > 0) %>%
dplyr::group_by(north) %>%
dplyr::mutate(season = factor(month.to.season(month, north),
levels = c('Spring', 'Summer', 'Fall', 'Winter'),
ordered = T)) %>%
ungroup()
}
get.data('precip.chirps.ww.monthly', overwrite = F,
load.precip.monthly,
precip.chirps.ww.monthly.path,
slide.all.precip)
precip.chirps.ww.seasonal.median <- precip.chirps.ww.monthly %>%
group_by(OBJECTID, year, season) %>%
summarize(precip.mm = sum(precip.mm)) %>%
ungroup() %>%
summarize(median = median(precip.mm)) %>%
pull(median)
precip.chirps.ww.monthly
precip.chirps.ww.event %>%
# Filter out sites with no precipitation data at all
filter(days_before==-1, window==7, !is.na(pctl)) %>%
# Filter out sites already excluded by landslide criteria
inner_join(slide.all.precip, by = 'OBJECTID') %>%
# Count numbers of zero and non-zero precipitation
group_by(precip = ifelse(pctl > 0, 'non.zero.precip', 'zero.precip')) %>%
tally() %>%
spread(key = precip, value = n) %>%
mutate(total = non.zero.precip + zero.precip,
percent.zero.precip = zero.precip / total)
filter.slide.precip <- function (slide.all.precip, precip.event) {
slide.all.precip %>%
dplyr::inner_join(precip.event %>%
filter(days_before==-1, window==7, !is.na(pctl)) %>%
dplyr::filter(days_before==-1, window==7, pctl > 0)) %>%
dplyr::select(OBJECTID, latitude, longitude, event_date, event_yday, season,
location_accuracy,
landslide_category, landslide_trigger, landslide_size,
landslide_setting, debris_flow)
}
get.data('slide', overwrite = F,
filter.slide.precip,
slide.all.precip,
precip.chirps.ww.event)
slide
filter.slide.ca <- function (slide.all.precip) {
slide.all.precip %>%
filter(is.in.ca.study(longitude, latitude)) %>%
mutate(region = as.factor(case_when(
latitude < 35.0 ~ 'South',
longitude < -121.25 ~ 'North',
TRUE ~ 'East')))
}
get.data('slide.ca', overwrite = F,
filter.slide.ca,
slide.all.precip)
slide %>% summary
## OBJECTID latitude longitude event_date
## Min. :388907 Min. :-46.77 Min. :-179.865 Min. :2007-01-05
## 1st Qu.:391763 1st Qu.: 10.33 1st Qu.: -96.159 1st Qu.:2010-11-30
## Median :394578 Median : 28.55 Median : -8.624 Median :2013-06-16
## Mean :394614 Mean : 23.09 Mean : -2.247 Mean :2013-05-15
## 3rd Qu.:397488 3rd Qu.: 38.95 3rd Qu.: 92.212 3rd Qu.:2016-03-10
## Max. :400281 Max. : 49.98 Max. : 179.402 Max. :2018-08-07
##
## event_yday season location_accuracy landslide_category
## Min. : 1.0 Spring:1310 exact: 423 landslide :3492
## 1st Qu.: 91.0 Summer:1798 1km :1573 mudslide :1326
## Median :172.0 Fall :1104 5km :2259 rock_fall : 226
## Mean :171.6 Winter:1101 10km :1058 complex : 103
## 3rd Qu.:244.0 debris_flow: 98
## Max. :366.0 (Other) : 64
## NA's : 4
## landslide_trigger landslide_size landslide_setting
## rain :1656 small :1735 above_road :1257
## continuous_rain: 561 medium :3160 natural_slope: 265
## downpour :3054 large : 326 urban : 178
## flooding : 42 very_large : 36 below_road : 121
## catastrophic: 3 above_river : 65
## NA's : 53 (Other) : 187
## NA's :3240
## debris_flow
## Mode :logical
## FALSE:3987
## TRUE :1326
##
##
##
##
filter.precip.ca <- function (precip.chirps.ww.event, slide.ca) {
precip.chirps.ca <- precip.chirps.ww.event%>%
dplyr::filter(OBJECTID %in% slide.ca$OBJECTID)
}
get.data('precip.chirps.ca.event', filter.precip.ca, overwrite = F,
precip.chirps.ww.event, slide.ca)
get.data('precip.chirps.ca.all', filter.precip.ca, overwrite = F,
precip.chirps.ww.all, slide.ca)
precip.chirps.ca.event %>%
# Filter out sites with no precipitation data at all
filter(days_before==-1, window==7, !is.na(pctl)) %>%
# Filter out sites already excluded by landslide criteria
inner_join(slide.all.precip, by = 'OBJECTID') %>%
# Count numbers of zero and non-zero precipitation
group_by(precip = ifelse(pctl > 0, 'non.zero.precip', 'zero.precip')) %>%
tally() %>%
spread(key = precip, value = n) %>%
mutate(total = non.zero.precip + zero.precip,
percent.zero.precip = zero.precip / total * 100)
load.precip.daymet.ca <- function (precip.daymet.ca.path) {
read_csv(precip.daymet.ca.path,
col_types = cols(X1 = col_datetime(format = "%Y-%m-%d %H:%M:%S"))) %>%
rename(date=X1) %>%
date.filter.precip %>%
tidy.precip
}
get.data('precip.daymet.ca.all', overwrite = F,
load.precip.daymet.ca,
precip.daymet.ca.path)
precip.daymet.ca.all
load.modis <- function (modis.paths, burn.date.min, burn.date.max) {
modis.date.raw <- dplyr::bind_rows(map(modis.paths, read_csv))
burn.record.length <- time_length(burn.date.max - burn.date.min, 'years')
modis.date.raw %>%
dplyr::right_join(slide %>% dplyr::select(OBJECTID, event_date, latitude),
by='OBJECTID') %>%
# Convert MODIS Burned Area value to date
dplyr::mutate(burn_date = make_date(year=floor(month / 1000)) + burn_date) %>%
dplyr::select(-month) %>%
# Separate out indivicual fires
dplyr::group_by(OBJECTID) %>%
dplyr::mutate(gap = as.integer(burn_date - lag(burn_date, default=min(burn_date))),
new.fire = gap > 30,
fire.id = cumsum(new.fire)) %>%
drop_na(fire.id) %>%
# Calculate burn fraction, fire ending date, and season for each fire
dplyr::group_by(OBJECTID, event_date, fire.id) %>%
summarize(fraction = sum(fraction),
burn_date = max(burn_date),
burn_season = yday.to.season(yday(burn_date), latitude)) %>%
# Calculate time relative to landslide
dplyr::mutate(burn_days_before = event_date - burn_date,
burn_years_before = time_length(burn_days_before, 'years')) %>%
# Compute fire information
dplyr::group_by(OBJECTID) %>%
dplyr::mutate(pre.slide = burn_date==max(burn_date[burn_days_before>=0]),
n = n(),
period = burn.record.length / n,
start_days_before = min(event_date) - burn.date.min,
start_years_before = time_length(start_days_before, 'years'),
end_days_before = max(event_date) - burn.date.max,
end_years_before = time_length(end_days_before, 'years')) %>%
dplyr::select(-event_date) %>%
dplyr::ungroup()
}
get.data('modis.date', overwrite = F,
load.modis,
modis.paths, modis.start.date, modis.end.date)
modis.date
classify.burned <- function (modis.date, slide) {
modis.date %>%
dplyr::filter(burn_years_before < 3, pre.slide) %>%
dplyr::right_join(slide, by='OBJECTID') %>%
dplyr::group_by(OBJECTID) %>%
summarize(burn.cumulative = sum(fraction)) %>%
dplyr::mutate(burn.cumulative = ifelse(is.na(burn.cumulative),
0,
burn.cumulative),
burned = ifelse(burn.cumulative > 0, 'burned', 'unburned')) %>%
dplyr::select(OBJECTID, burned, burn.cumulative) %>%
dplyr::ungroup()
}
get.data('modis', overwrite = F,
classify.burned,
modis.date,
slide)
modis %>%
group_by(burned) %>%
tally() %>%
mutate(percent = n / nrow(slide) * 100)
cluster.global <- function (slide, agclusts, k=30) {
# Select clusters with more than 100 landslides, and label
slide %>%
dplyr::select(OBJECTID, latitude, longitude) %>%
dplyr::mutate(region = agclusts %>% cutree(k)) %>%
dplyr::group_by(region) %>%
mutate(region = as_factor(region)) %>%
dplyr::filter(n() > 100) %>%
dplyr::ungroup() %>%
dplyr::mutate(region = fct_recode(as_factor(region),
'Central Europe' = '1',
'Malaysia' = '2',
'Western US/Canada' = '3',
'East Brazil' = '4',
'Southeast Asia' = '5',
'South India' = '6',
'Eastern US' = '8',
'Himalayas' = '9',
'Southeast Asia' = '10',
'Caribbean/Venezuela' = '11',
'Central America' = '14'))
}
# Cluster landslides based on location using AGNES algorithm
get.data('agclusts',
cluster::agnes,
slide %>%
dplyr::select(latitude, longitude))
get.data('clusters.global', overwrite = F,
cluster.global,
slide.all.precip, agclusts)
clusters.global %>%
ggplot(aes(longitude, latitude)) +
theme_global_map +
geom_point(aes(color = region), size=1) +
scale_color_brewer('Region', palette='Paired') +
guides(color=guide_legend(override.aes = list(size=5)))
split.regions.us <- function (agclusts, slide, k=6) {
# Crop cluster tree, select clusters with more than 100 landslides, and label
slide %>%
dplyr::select(OBJECTID, latitude, longitude) %>%
dplyr::mutate(region = agclusts %>% cutree(k)) %>%
mutate(region = as_factor(region)) %>%
dplyr::mutate(region = fct_recode(as_factor(region),
'Intermountain West' = '3',
'Intermountain West' = '4',
'Intermountain West' = '5',
'Pacific Northwest' = '1',
'California' = '2')) %>%
dplyr::group_by(region) %>%
dplyr::filter(n() > 100) %>%
dplyr::ungroup() %>%
droplevels()
}
get.data('agclusts.us', overwrite = F,
function (coord) cluster::agnes(coord, stand=T),
clusters.global %>%
dplyr::filter(region=='Western US/Canada') %>%
dplyr::select(latitude, longitude))
get.data('clusters.us', overwrite = F,
split.regions.us,
agclusts.us,
clusters.global %>%
dplyr::filter(region=='Western US/Canada'))
clusters.us %>%
ggplot(aes(longitude, latitude)) +
geom_point(aes(color = region), size=1) +
scale_color_brewer(palette='Set2') +
guides(color=guide_legend(override.aes = list(size=5))) +
theme_map
merge.regions <- function (clusters.global, clusters.us) {
# Merge US and global regions
clusters.global %>%
left_join(clusters.us %>% dplyr::select(OBJECTID, region),
by='OBJECTID') %>%
dplyr::mutate(region = as_factor(ifelse(!is.na(region.y),
as.character(region.y),
as.character(region.x)))) %>%
dplyr::select(-starts_with('region.')) %>%
filter(region != 'Western US/Canada') %>%
droplevels()
}
get.data('regions.all', overwrite = F,
merge.regions,
clusters.global, clusters.us)
get.regions.info <- function (regions.all, slide) {
regions.all %>%
bind_rows(slide %>% dplyr::mutate(region = 'All landslides')) %>%
dplyr::inner_join(modis, by = 'OBJECTID') %>%
group_by(region) %>%
dplyr::group_by(region, burned) %>%
dplyr::summarize(n = n()) %>%
tidyr::spread(key = burned, value = n) %>%
dplyr::mutate(n = burned + unburned,
percent.burned = burned / n * 100) %>%
ungroup
}
get.data('regions.info', overwrite = F,
get.regions.info,
regions.all, slide)
regions.info
filter.and.merge.regions <- function (regions.all, regions.info, slide) {
# Combine nearby regions
regions.all %>%
filter(OBJECTID %in% slide$OBJECTID) %>%
bind_rows(slide %>% dplyr::mutate(region = 'All landslides')) %>%
dplyr::left_join(regions.info) %>%
mutate(region = fct_recode(region,
'Central America' = 'Caribbean/Venezuela')) %>%
dplyr::filter(region != 'East Brazil',
percent.burned > 4) %>%
droplevels() %>%
dplyr::mutate(region = factor(region, region.plot.order, ordered = T)) %>%
dplyr::select(OBJECTID, region)
}
get.data('slide.regions', overwrite = F,
filter.and.merge.regions,
regions.all,
regions.info,
slide)
regions.tally <- slide.regions %>%
group_by(region) %>%
tally
regions.tally
A panelled maps showing landslide regions (a and b), whether or not each landslide took place on a burned site (c and d), and the precipitation percentile of the landslide-triggering storm (e and f). Inset plots show, for each region, the distribution of burned fractions of the landslide sites (a and b), location accuracy of landslides split by fire history (c and d), and the seasonal total precipitation climatology (e and f). ### Regions map
gg.region <- slide.regions %>%
left_join(slide, by='OBJECTID') %>%
mutate(region = fct_relevel(region, 'All landslides', after = Inf)) %>%
arrange(desc(region)) %>%
ggplot(aes(longitude, latitude)) +
theme_global_map +
geom_point(aes(color = region), size=0.02) +
labs(color = '') +
scale_color_region +
guides(color = guide_legend(override.aes = list(size = 2))) +
theme_panel_legend
gg.region
gg.burned <- slide %>%
dplyr::left_join(modis) %>%
dplyr::arrange(desc(burned)) %>%
ggplot() +
theme_global_map +
#geom_point(data = slide.regions %>%
# dplyr::filter(region !='All landslides') %>%
# dplyr::left_join(slide, by='OBJECTID'),
# aes(x = longitude, y = latitude, fill = region),
# size = 2, alpha = .01, shape = 21, stroke = NA) +
geom_point(aes(x=longitude, y=latitude, color=burned),
size=.02) +
scale_color_burned +
#scale_fill_region +
guides(color = guide_legend(override.aes = list(size = 2),
title = NULL),
fill = guide_none()) +
theme_panel_legend
gg.burned
gg.precip <- slide %>%
left_join(precip.chirps.ww.event %>%
filter(days_before==-1, window==7),
by='OBJECTID') %>%
filter(pctl > 0) %>%
arrange(desc(pctl)) %>%
ggplot() +
theme_global_map +
geom_point(aes(x=longitude, y=latitude, color=pctl),
size=.02) +
scale_color_precip +
theme_panel_legend
gg.precip
min.burn.cumulative <- modis %>%
filter(burn.cumulative > 0) %>%
pull(burn.cumulative) %>%
min
max.burn.cumulative <- modis %>%
pull(burn.cumulative) %>%
max
plot.burned.fraction <- function (reg, slide.regions, modis) {
slide.regions %>%
filter(region==reg) %>%
left_join(modis, by = 'OBJECTID') %>%
filter(burn.cumulative > 0) %>%
ggplot() +
geom_density(aes(x = burn.cumulative), fill = 'black') +
scale_x_log10() +
labs(x = 'Burned fraction', y = '') +
scale_color_region +
scale_fill_region +
coord_cartesian(expand=F, xlim = c(min.burn.cumulative, max.burn.cumulative)) +
theme_bw() +
theme(legend.position = 'none',
axis.text.x = element_text(angle = 90, vjust = 0.5, size = 4),
axis.title = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 5, margin=margin(0,0,0,0))) +
ggtitle(str_c(reg, '\nn=', pull(filter(regions.tally, region==reg))))
}
gg.burned.fraction <- region.plot.order %>%
set_names %>%
map(plot.burned.fraction, slide.regions, modis)
cowplot::plot_grid(plotlist = gg.burned.fraction, nrow = 1)
plot.accuracy <- function (reg, slide.regions, slide, modis) {
slide.regions %>%
left_join(slide, by='OBJECTID') %>%
left_join(modis, by='OBJECTID') %>%
filter(region==reg) %>%
mutate(location_accuracy = fct_recode(location_accuracy, '0km' = 'exact')) %>%
drop_na(location_accuracy) %>%
ggplot() +
geom_bar(aes(x = location_accuracy, fill=burned,
group = interaction(region, burned),
y = ..prop..),
position=position_dodge2(preserve = 'single')) +
scale_fill_burned +
coord_cartesian(expand=F, ylim = c(0, 0.6)) +
theme_bw() +
labs(y = '', x = '', fill = 'Location\naccuracy') +
theme(legend.position = 'none',
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1,
size = 4),
axis.title = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 5, margin=margin(0,0,0,0))) +
ggtitle(str_c(reg, '\nn=', pull(filter(regions.tally, region==reg))))
}
gg.accuracy <- region.plot.order %>%
set_names %>%
map(plot.accuracy, slide.regions, slide, modis)
cowplot::plot_grid(plotlist = gg.accuracy, nrow = 1)
plot.climatology <- function (reg, precip.monthly, slide.regions, med) {
precip.seasonal <- precip.monthly %>%
group_by(OBJECTID, year, season) %>%
summarize(precip.mm = sum(precip.mm))
ymin <- 1
ymax <- max(precip.seasonal$precip.mm, na.rm=T)
precip.seasonal %>%
dplyr::right_join(slide.regions %>% filter(region==reg)) %>%
drop_na(season) %>%
ggplot() +
geom_violin(aes(y = precip.mm, x = season, group = season),
fill='black', color=NA, adjust = 2) +
geom_hline(aes(yintercept = med), color = 'grey', size = .3)+
scale_y_log10(breaks = c(1, 10, 100, 1000)) +
scale_x_discrete(labels = c('Spring' = 'MAM', 'Summer' = 'JJA',
'Fall' = 'SON', 'Winter' = 'DJF'),
expand = ggplot2::expand_scale(add = c(1, 0))) +
coord_cartesian(expand=F, ylim = c(ymin, ymax)) +
theme_bw() +
theme(legend.position = 'none',
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 4,
margin=margin(0,0,0,0)),
axis.text.y = element_text(size = 4, hjust = 1, vjust = 0.5),
axis.title = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 5, margin=margin(0,0,0,0))) +
ggtitle(str_c(reg, '\nn=', pull(filter(regions.tally, region==reg))))
}
gg.climatology <- region.plot.order %>%
purrr::set_names() %>%
purrr::map(plot.climatology,
precip.chirps.ww.monthly,
slide.regions,
precip.chirps.ww.seasonal.median)
cowplot::plot_grid(plotlist = gg.climatology, nrow = 1)
coord_america <- coord_equal(xlim = c(-140, -50), ylim = c(-8, 53), expand = F)
coord_asia <- coord_equal(xlim = c(50, 140), ylim = c(-12, 49), expand = F)
# Clip maps
gg.america.region <- gg.region + theme_panel + coord_america
gg.asia.region <- gg.region + theme_panel + coord_asia
gg.america.burned <- gg.burned + theme_panel + coord_america
gg.asia.burned <- gg.burned + theme_panel + coord_asia
gg.america.precip <- gg.precip + theme_panel + coord_america
gg.asia.precip <- gg.precip + theme_panel + coord_asia
# Pull legends
legend.regions <- cowplot::get_legend(gg.region)
legend.burned <- cowplot::get_legend(gg.burned)
legend.precip <- cowplot::get_legend(gg.precip)
# Region midpoints for arrows
region.centroids <- slide.regions %>%
left_join(slide) %>%
group_by(region) %>%
summarize(y = mean(latitude), x = mean(longitude))
# Add insets
inset.w <- .2
inset.h <- .41
inset.arrow <- grid::arrow(length = unit(0.03, "npc"), type = 'open')
# Regions
gg.america.region.inset <- cowplot::ggdraw(gg.america.region) +
# Arrows
cowplot::draw_grob(grid::grid.lines(x = c(.1, .25), y = c(.8, .8),
arrow=inset.arrow)) +
cowplot::draw_grob(grid::grid.lines(x = c(.1, .275), y = c(.4, .65),
arrow=inset.arrow)) +
cowplot::draw_grob(grid::grid.lines(x = c(.4, .4), y = c(.375, .575),
arrow=inset.arrow)) +
cowplot::draw_grob(grid::grid.lines(x = c(.825, .675), y = c(.675, .3),
arrow=inset.arrow)) +
# Inset plots
cowplot::draw_plot(gg.accuracy$`Pacific Northwest`,
x = 0.01, y = 1,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.accuracy$California,
x = 0.01, y = .6,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.accuracy$`Intermountain West`,
x = .275, y = .4,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.accuracy$`Central America`,
x = .775, y = .8,
width = inset.w, height = inset.h, vjust = 1)
gg.asia.region.inset <- cowplot::ggdraw(gg.asia.region) +
#Arrows
cowplot::draw_grob(grid::grid.lines(x = c(.3, .4), y = c(.3, .6),
arrow=inset.arrow)) +
cowplot::draw_grob(grid::grid.lines(x = c(.825, .8), y = c(.55, .4),
arrow=inset.arrow)) +
# Inset plots
cowplot::draw_plot(gg.accuracy$Himalayas,
x = .15, y = 0.4,
width = inset.w, height= inset.h, vjust = 1) +
cowplot::draw_plot(gg.accuracy$`Southeast Asia`,
x = .775, y = .935,
width = inset.w, height = inset.h, vjust = 1)
# Fire
gg.america.burned.inset <- cowplot::ggdraw(gg.america.burned) +
cowplot::draw_plot(gg.burned.fraction$`Pacific Northwest`,
x = 0, y = 1,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.burned.fraction$California,
x = 0, y = .6,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.burned.fraction$`Intermountain West`,
x = .275, y = .4,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.burned.fraction$`Central America`,
x = .775, y = .8,
width = inset.w, height = inset.h, vjust = 1)
gg.asia.burned.inset <- cowplot::ggdraw(gg.asia.burned) +
cowplot::draw_plot(gg.burned.fraction$Himalayas,
x = .15, y = 0.4,
width = inset.w, height= inset.h, vjust = 1) +
cowplot::draw_plot(gg.burned.fraction$`Southeast Asia`,
x = .775, y = .935,
width = inset.w, height = inset.h, vjust = 1)
# Precipitation
inset.w.p = .23
inset.h.p = .41
gg.america.precip.inset <- cowplot::ggdraw(gg.america.precip) +
cowplot::draw_plot(gg.climatology$`Pacific Northwest`,
x = 0, y = 1.01,
width = inset.w.p, height = inset.h.p, vjust = 1) +
cowplot::draw_plot(gg.climatology$California,
x = .085, y = .64,
width = inset.w.p, height = inset.h.p, vjust = 1) +
cowplot::draw_plot(gg.climatology$`Intermountain West`,
x = .245, y = .395,
width = inset.w.p, height = inset.h.p, vjust = 1) +
cowplot::draw_plot(gg.climatology$`Central America`,
x = .71, y = .8,
width = inset.w.p, height = inset.h.p, vjust = 1)
gg.asia.precip.inset <- cowplot::ggdraw(gg.asia.precip) +
cowplot::draw_plot(gg.climatology$Himalayas,
x = .15, y = 0.4,
width = inset.w.p, height= inset.h.p, vjust = 1) +
cowplot::draw_plot(gg.climatology$`Southeast Asia`,
x = .7, y = .935,
width = inset.w.p, height = inset.h.p, vjust = 1)
panel.labels <- tibble(
label = as_factor(c('Location accuracy distribution\nfor burned and unburned sites by region',
'Locations of burned sites and\nregional distribution of burned fraction',
str_c('Landslide-triggering precipitation percentile\nand',
'regional seasonal precipitation (mm)')))
)
gg.panel.labels <- ggplot(panel.labels) +
facet_wrap(~label, ncol = 1, switch = T) +
theme(panel.background = element_blank(),
strip.text = element_text(face="bold", size=7))
gg.panel.maps <- cowplot::plot_grid(
cowplot::plot_grid(gg.america.region.inset, gg.asia.region.inset,
labels = c('a', 'b'), label_x = 1, hjust = 1, vjust = 1),
legend.regions,
cowplot::plot_grid(gg.america.burned.inset, gg.asia.burned.inset,
labels = c('c', 'd'), label_x = 1, hjust = 1, vjust = 1),
legend.burned,
cowplot::plot_grid(gg.america.precip.inset, gg.asia.precip.inset,
labels = c('e', 'f'), label_x = 1, hjust = 1, vjust = 1),
legend.precip,
ncol = 1, rel_heights = c(1, .3, 1, .2, 1, .3))
gg.panel <- cowplot::plot_grid(gg.panel.labels, gg.panel.maps,
nrow = 1, rel_widths = c(.5, 6.5))
gg.panel
ggsave(filename = "../figures/panel_map.png", width=7, height=8)
This figure shows the timing of fires relative to landslides for all events that took place at a burned site.
calc.fire.timing <- function (modis.date, modis, slide) {
modis.date %>%
dplyr::filter(burn_days_before >= 0) %>%
dplyr::group_by(OBJECTID) %>%
dplyr::filter(burn_days_before == min(burn_days_before)) %>%
right_join(modis %>% dplyr::filter(burned=='burned'), by='OBJECTID') %>%
inner_join(slide.regions, by='OBJECTID') %>%
left_join(slide %>% dplyr::select(OBJECTID, event_yday, season), by='OBJECTID') %>%
dplyr::mutate(burn_yday = yday(burn_date),
burn_season = factor(burn_season,
levels = c('Spring', 'Summer', 'Fall', 'Winter'),
ordered = T))
}
get.data('fire.timing', overwrite = F,
calc.fire.timing,
modis.date, modis, slide)
Most landslides take place within one year of the fire. California stands out as having over a quarter of the landslides taking place more than a year after the fire.
fire.timing %>%
dplyr::mutate(gt.1year = ifelse(burn_days_before > 2*365,
'gt.1year', 'lt.1year')) %>%
dplyr::group_by(region, gt.1year) %>%
dplyr::tally() %>%
tidyr::spread(gt.1year, n, drop=F) %>%
dplyr::mutate(prop.gt.1year = gt.1year / (gt.1year + lt.1year))
fire.to.landslide.x.breaks <- seq(season.breaks[4] - 4 * 365, season.breaks[3],
365/4)
fire.to.landslide.x.labels <- map(c(str_c('-', c(3, 2, 1), 'y '), ''),
~c(str_c(., 'Winter'),
'Spring', 'Summer', 'Fall')) %>% as_vector
fire.facet.labels <- fire.timing %>%
dplyr::group_by(region) %>%
dplyr::mutate(wrap_event_yday = ifelse(event_yday > season.breaks['winter'],
event_yday - 365, event_yday)) %>%
dplyr::ungroup() %>%
dplyr::mutate(ymin = min(wrap_event_yday - burn_days_before),
ymax = max(wrap_event_yday)) %>%
dplyr::group_by(region) %>%
dplyr::summarize(y = n() * 0.8, x = min(ymin) + (max(ymax) - min(ymin)) * 0.1) %>%
dplyr::mutate(label = c('a', 'b', 'c', 'd', 'e', 'f', 'g'))
fire.timing %>%
dplyr::group_by(region) %>%
dplyr::mutate(wrap_event_yday = ifelse(event_yday > season.breaks['winter'],
event_yday - 365, event_yday)) %>%
dplyr::arrange(wrap_event_yday - burn_days_before) %>%
dplyr::mutate(i = row_number()) %>%
ggplot() +
facet_wrap(vars(region), scales = 'free_y', ncol = 2, shrink=F) +
geom_segment(aes(x = wrap_event_yday - burn_days_before, xend = wrap_event_yday,
y = i, yend = i,
color = burn_season)) +
geom_point(aes(x = wrap_event_yday, y = i, shape='Landslide'), size=0.5) +
geom_vline(data = tibble(x=seq(334-365*3, 365, 365)),
aes(xintercept = x, size = 'start_of_winter')) +
geom_vline(aes(xintercept = 334 - 365, size = 'year_of')) +
geom_rug(aes(y = i, color = burn_season)) +
geom_rug(aes(x = event_yday - burn_days_before), sides = 't', length = unit(0.05, "npc")) +
geom_label(data = fire.facet.labels,
aes(x = x, y = y, label = label)) +
theme_bw() +
scale_x_continuous('Landslides and preceding fires relative to the year of the landslide',
minor_breaks = NULL,
breaks = as_vector(fire.to.landslide.x.breaks),
labels = fire.to.landslide.x.labels) +
scale_color_manual('Fire season:',
values = c('Spring'='#19E667', 'Summer'='#E61998',
'Fall'='#E6CD19', 'Winter'='#1932E6'),
guide = guide_legend(nrow=1, order = 2)) +
scale_y_continuous('Post-wildfire landslides ordered by fire timing relative to landslide',
expand = c(0.1, 0.1)) +
scale_size_manual('Relative timing:',
values = c('year_of' = 1, 'start_of_winter' = 0.5),
labels = c('year_of' = 'Start of year of landslide',
'start_of_winter' = 'Start of winter (December 1)'),
guide = guide_legend(nrow=2, order = 3)) +
guides(shape = guide_legend(title = NULL, override.aes = c(size = 1.5), order = 1)) +
theme(panel.grid.major.x = element_line(color = 'grey', size = 0.5),
axis.text.x = element_text(hjust = 1, vjust = 1, angle = 90),
axis.ticks = element_blank(),
axis.text.y = element_blank(),
legend.position = c(.75, .04),
legend.direction = 'vertical',
legend.box = 'vertical',
legend.margin = unit(5, 'pt'),
legend.box.margin = margin(5, 5, 5, 5),
legend.box.background = element_rect(color = 'black', size = 2),
legend.key.height = unit(10, 'pt')
)
ggsave('../figures/fire_to_landslide.png', width = 6.5, height = 8)
To validate the exclusion of zero-precipitation landslides, a comparison of CHIRPS and Daymet in California and Nevada only. Many of the zero-precipitation landslides measured by CHIRPS are not the same as the ones measured by Daymet, suggesting that these are measurement errors.
Approximately 14% of the landslides in California and Nevada were excluded on this basis.
join.chirps.daymet <- function (precip.chirps, precip.daymet) {
precip.chirps %>%
dplyr::inner_join(precip.daymet %>% dplyr::select(-days_before, -yday),
by=c('OBJECTID', 'date', 'window'),
suffix=c('.chirps', '.daymet'))
}
get.data('precip.chirps.daymet.ca.event', join.chirps.daymet, overwrite = F,
precip.chirps.ca.event,
precip.daymet.ca.all)
get.data('precip.chirps.daymet.ca.all', join.chirps.daymet, overwrite = F,
precip.chirps.ca.all, precip.daymet.ca.all)
precip.chirps.ca.event %>%
filter(days_before==-1, window==7, !is.na(pctl)) %>%
group_by(precip = ifelse(pctl > 0, 'non.zero.precip', 'zero.precip')) %>%
tally() %>%
spread(key = precip, value = n) %>%
mutate(total = non.zero.precip + zero.precip,
percent.zero.precip = zero.precip / total)
precip.chirps.daymet.ca.all %>%
tidyr::drop_na(pctl.chirps, pctl.daymet) %>%
ggplot() +
geom_point(aes(x=pctl.chirps, y=pctl.daymet, color = 'all'),
alpha=0.2,shape = 16) +
geom_point(data = precip.chirps.daymet.ca.event %>%
dplyr::filter(window == 7, days_before == -1 ) %>%
dplyr::mutate(to.remove = ifelse(pctl.chirps==0, 'exclude', 'keep')),
aes(x=pctl.chirps, y=pctl.daymet, color=to.remove)) +
coord_equal() +
theme_bw() +
scale_color_manual('',
values = c('keep'='black', 'exclude'='blue', 'all'='gray'),
labels = c('keep'='Included\nday-of-landslide\npreciptiation',
'exclude'='Excluded\nday-of-landlside\nprecipitation',
'all'='Precipitation\nfrom all dates')) +
guides(color = guide_legend(keyheight = grid::unit(4, 'char'))) +
labs(x='CHIRPS 7-day Precipitation Percentile',
y='Daymet 7-day Precipitation Percentile')
ggsave('../figures/chirps_vs_daymet.png', width=6, height=4, dpi=600)
Comparison of antecedent precipitation time series among burned and unburned sites by region. Mann-Whitney tests are used to establish whether or not differences are statistically significant.
calc.burn.wilcox <- function (precip, pre=30, post=10, alternative='two.sided',
source=NA, groupby=c('burned', 'days_before')) {
precip %>%
dplyr::filter(days_before <= pre, days_before >= -post) %>%
dplyr::group_by_at(groupby) %>%
calc.wilcox(alternative, source)
}
calc.wilcox <- function (precip, alternative='two.sided', source=NA) {
wilcox <- precip %>%
tidyr::nest() %>%
spread(key = 'burned', value = 'data') %>%
dplyr::mutate(wilcox = map2(burned, unburned,
~wilcox.test(.x$pctl, .y$pctl,
alternative=alternative))) %>%
dplyr::mutate(wilcox.glance = map(wilcox, glance)) %>%
tidyr::unnest(wilcox.glance) %>%
dplyr::mutate(signif.95 = p.value <= 0.05) %>%
gather(key = 'burned', value = 'data', burned, unburned) %>%
tidyr::unnest(data)
if (!is.na(source)) wilcox <- wilcox %>% dplyr::mutate(label=source)
wilcox
}
magnitude.regional.data <- function (slide.regions, precip,
pre = 30, post = 10, w = 7) {
precip %>%
dplyr::filter(days_before < pre, days_before > -post, window==w) %>%
right_join(slide.regions, by = 'OBJECTID') %>%
dplyr::left_join(modis, by = 'OBJECTID') %>%
dplyr::filter(pctl > 0) %>%
droplevels() %>%
dplyr::select(region, burned, days_before, pctl) %>%
calc.burn.wilcox(pre = pre, post = post, alternative = 'less',
groupby = c('region', 'burned', 'days_before'))
}
get.data('burn.wilcox.regions.event', overwrite = F,
magnitude.regional.data,
slide.regions,
precip.chirps.ww.event)
facet.levels <- c('Rolling cumulative\nprecipitation percentile',
'Mann-Whitney p-value')
plot.burn.wilcox <- function (burn.wilcox.regions.event, slide.regions) {
region.tally <- slide.regions %>%
left_join(modis) %>%
dplyr::group_by(burned=burn.cumulative > 0, region) %>%
tally
burn.wilcox.regions.event %>%
dplyr::mutate(signif.95 = ifelse(signif.95, 'significant', 'not.significant')) %>%
ggplot() +
facet_wrap(region~., ncol=1, strip.position = 'right') +
geom_boxplot(aes(x = days_before,
group=interaction(days_before, burned),
y = pctl,
fill = burned,
alpha = interaction(signif.95, burned)),
outlier.shape=21) +
annotate("rect", xmin = -0.5, xmax = 0.5, ymin = 0, ymax = 1, alpha = .4) +
geom_label(data=region.tally,
aes(x=-15, y=0.1+0.25*burned,
label=str_c(ifelse(burned, 'Burned: ', 'Unburned: '), n)),
hjust='right', label.size = 0) +
geom_label(data = tibble(region = region.plot.order,
label = c('a', 'b', 'c', 'd', 'e', 'f', 'g')) %>%
dplyr::mutate(region = factor(region, region.plot.order, ordered = T)),
aes(x = -15, y = 0.85, label = label), hjust = 'right') +
scale_x_reverse(limits = c(31, -16), expand = expansion()) +
scale_fill_manual('',
values = c('burned' = '#d95f02', unburned = '#7570b3'),
labels = c('burned'='Burned', 'unburned'='Unburned'),
guide = guide_legend(direction = 'horizontal')) +
scale_alpha_manual('Statistical significance at\n95% confidence',
values=c('not.significant.burned' = 0.4,
'not.significant.unburned' = 0.4,
'significant.burned' = 1,
'significant.unburned' = 1),
labels=c('not.significant.burned' = '',
'not.significant.unburned' = 'Not Significant',
'significant.burned' = '',
'significant.unburned' = 'Significant'),
breaks=c('not.significant.burned',
'not.significant.unburned',
'significant.burned',
'significant.unburned'),
guide=guide_legend(override.aes=list(fill=c('not.significant.burned' = '#d95f02',
'not.significant.unburned' = '#7570b3',
'significant.burned' = '#d95f02',
'significant.unburned' = '#7570b3')),
direction = 'horizontal')) +
theme_bw() +
labs(x='Days before the event', y='Rolling cumulative precipitation percentile') +
theme(legend.position = 'bottom', legend.direction = 'vertical')
}
burn.wilcox.regions.event %>%
plot.burn.wilcox(slide.regions)
ggsave('../figures/magnitude_regional.png', width = 11, height = 11)
Bootstrap samples are used here to develop baseline precipitation values against which to compare landslide-triggering precipitation. This method is an improvement over the original Wilcox tests because sample-size disparities between regions are reduced.
sample.bootstrap <- function(precip.subset, precip.all, n, samples, burned, region) {
precip.all %>%
dplyr::filter(burned==burned, region==region) %>%
dplyr::mutate(precip.yday = yday(date)) %>%
# dplyr::select precipitation from same location and date of year as subset
dplyr::inner_join(precip.subset %>%
dplyr::select(OBJECTID, event_yday),
by=c('OBJECTID'='OBJECTID',
'precip.yday'='event_yday')) %>%
# Sample precipitation matching the site/day in bulk
dplyr::group_by(OBJECTID) %>%
tidyr::drop_na(pctl) %>%
dplyr::sample_n(n, replace = T) %>%
# Randomly split into individual samples
dplyr::mutate(i = sample(1:samples, n, replace = T ))
}
calc.burn.bootstrap <- function (precip.event, precip.all,
samples=100, per.subset=10000,
pre=30, post=10, w=7) {
groupby = c('region', 'burned', 'days_before')
precip.all <- precip.all %>%
dplyr::filter(window==w) %>%
dplyr::select(-window)
precip.event %>%
rename(precip.date = date, precip.yday = yday) %>%
dplyr::filter(days_before <= pre, days_before >= -post, window==w) %>%
dplyr::select_at(c(groupby, 'OBJECTID', 'event_yday')) %>%
# Tally each subset/site/day combination
dplyr::group_by_at(groupby) %>%
add_tally(name = 'sites.in.subset') %>%
# Determine number to draws from each site that will
# roughly match sample sizes among subsets
# Round up - number of draws must be an integer
dplyr::mutate(n.per.site = ceiling(per.subset / sites.in.subset)) %>%
# Take bootstrap samples for each subset/site/day combination
dplyr::group_by_at(c(groupby, 'n.per.site', 'sites.in.subset')) %>%
tidyr::nest() %>%
dplyr::mutate(bootstrap.samples = pmap(list(data, n.per.site,
burned, region),
~sample.bootstrap(..1, precip.all,
..2, samples,
..3, ..4)))
}
join.details <- function (precip, slide, slide.regions, modis) {
precip %>%
dplyr::inner_join(slide %>% dplyr::select(OBJECTID, event_yday),
by='OBJECTID') %>%
dplyr::right_join(slide.regions, by = 'OBJECTID') %>%
dplyr::left_join(modis %>% dplyr::select(OBJECTID, burned),
by='OBJECTID')
}
get.bootstrap <- function (reg='All landslides',
name = 'precip.chirps.bootstrap.all',
ovw=F) {
get.data(name, overwrite = ovw,
calc.burn.bootstrap,
precip.chirps.ww.event %>%
join.details(slide, slide.regions %>% filter(region == reg), modis),
precip.chirps.ww.all %>%
join.details(slide, slide.regions %>% filter(region == reg), modis))
}
get.bootstrap('California', 'precip.chirps.bootstrap.cal', ovw = F)
get.bootstrap('Intermountain West', 'precip.chirps.bootstrap.imw', ovw = F)
get.bootstrap('Pacific Northwest', 'precip.chirps.bootstrap.pnw', ovw = F)
get.bootstrap('Central America', 'precip.chirps.bootstrap.cam', ovw = F)
get.bootstrap('Himalayas', 'precip.chirps.bootstrap.him', ovw = F)
get.bootstrap('Southeast Asia', 'precip.chirps.bootstrap.sea', ovw = F)
get.bootstrap(ovw = F)
precip.chirps.bootstrap <- precip.chirps.bootstrap.all %>%
dplyr::bind_rows(precip.chirps.bootstrap.cal) %>%
dplyr::bind_rows(precip.chirps.bootstrap.imw) %>%
dplyr::bind_rows(precip.chirps.bootstrap.pnw) %>%
dplyr::bind_rows(precip.chirps.bootstrap.cam) %>%
dplyr::bind_rows(precip.chirps.bootstrap.him) %>%
dplyr::bind_rows(precip.chirps.bootstrap.sea)
precip.chirps.bootstrap %>%
dplyr::filter(days_before==0) %>%
dplyr::select(-data) %>%
dplyr::mutate(site.total = map_int(bootstrap.samples, ~nrow(ungroup(.))),
n.per.site.real = site.total / sites.in.subset,
avg.sample.size = map_dbl(bootstrap.samples,
~group_by(., i) %>%
summarize(n = n()) %>%
pull(n) %>% mean)) %>%
dplyr::select(-bootstrap.samples)
mann.whitney.p.value <- function (subset, bootstrap) {
bootstrap %>%
dplyr::group_by(i) %>%
tidyr::nest() %>%
# test if bootstrap samples are less than pre-landslide precipitation
dplyr::mutate(wilcox = map(data,
~wilcox.test(pull(., pctl),
pull(subset, pctl),
alternative='less'))) %>%
dplyr::mutate(wilcox.glance = map(wilcox, glance)) %>%
dplyr::select(-data, -wilcox) %>%
tidyr::unnest(wilcox.glance)
}
subset.precip <- function (precip, regions, modis, w=7, pre=30, post=10) {
precip %>%
dplyr::filter(window==7, days_before <= pre, days_before >= -post) %>%
dplyr::select(OBJECTID, pctl, days_before) %>%
dplyr::right_join(regions, by='OBJECTID') %>%
dplyr::inner_join(modis, by='OBJECTID') %>%
dplyr::group_by(region, burned, days_before) %>%
tidyr::nest()
}
get_subset <- function (reg, bur, day, subsets) {
row <- subsets %>%
dplyr::filter(region==as.character(reg),
burned==bur,
days_before==day)
row$data[[1]]
}
get.data('precip.chirps.ww.event.subsets',
subset.precip,
precip.chirps.ww.event,
slide.regions,
modis)
summarize.bootstrap <- function (bootstrap.df, subsets, objective.function) {
bootstrap.df %>%
dplyr::mutate(objective = purrr::pmap(list(region,
burned,
days_before,
bootstrap.samples),
~objective.function(get_subset(..1,
..2,
..3,
subsets),
..4))) %>%
dplyr::select(-bootstrap.samples) %>%
tidyr::unnest(objective) %>%
dplyr::ungroup()
}
get.data('precip.chirps.bootstrap.wilcox', overwrite = F,
summarize.bootstrap,
precip.chirps.bootstrap %>% dplyr::select(-data),
precip.chirps.ww.event.subsets,
mann.whitney.p.value)
precip.chirps.bootstrap.wilcox %>%
dplyr::select(region, burned, days_before, p.value)
The bootstrap plot also shows the day-of-landslide kernel density estimates for all bootstrap samples against the landslide-triggering precipitation to illustrate the method. ### Plot bootstrap sampling
plot.burn.bootstrap <- function (bootstrap, var, name) {
bootstrap %>%
ggplot() +
facet_wrap(.~region, ncol=1, strip.position='right') +
geom_vline(aes(xintercept=0)) +
geom_boxplot(aes(x = days_before,
group = interaction(days_before, burned),
y = {{ var }},
fill = burned,
color = burned),
outlier.shape=21,
alpha = 0.6,
position = position_dodge(width=0.4, preserve = 'total')) +
annotate("rect", xmin = -11, xmax = -0.1, ymin = 0.05, ymax = .9999,
fill = 'white', alpha = 0.4) +
annotate("rect", xmin = 0.1, xmax = 31, ymin = 0.05, ymax = .9999,
fill = 'white', alpha = 0.4) +
geom_hline(aes(yintercept=0.05), linetype='dashed') +
geom_label(data = tibble(region = region.plot.order,
label = c('a', 'b', 'c', 'd', 'e', 'f', 'g')) %>%
dplyr::mutate(region = factor(region, region.plot.order, ordered = T)),
aes(x = -11.5, y = 0.9999, label = label),
hjust = 'right', vjust = 'bottom') +
scale_x_reverse(expand = expansion(), limits = c(31, -12)) +
scale_y_continuous(trans = probit_trans(),
breaks = bootstrap.y.breaks,
labels = bootstrap.y.labels) +
scale_fill_burned +
scale_color_burned +
coord_trans(y = 'reverse',
ylim = c(10^-35, .9999)) +
theme_bw() +
labs(x='Days before the landslide', y=name) +
theme(legend.position = 'bottom',
panel.grid.minor.y = element_blank())
}
bootstrap.y.breaks <- c(10^-30, 10^-10, 0.05, 0.5, 0.9999)
bootstrap.y.labels <- c('1e-30', '1e-10', '0.05', '0.5', '1')
plot.bootstrap.wilcox <- precip.chirps.bootstrap.wilcox %>%
plot.burn.bootstrap(p.value,
'p-value of Wilcox tests between bootstrap samples and pre-landslide precipitation')
plot.bootstrap.wilcox
plot.bootstrap.sample.kde <- precip.chirps.bootstrap %>%
dplyr::filter(days_before==0) %>%
dplyr::select(-data) %>%
tidyr::unnest() %>%
ggplot() +
facet_grid(cols = vars(burned), rows = vars(region), scales = 'free_y') +
geom_density(aes(x = pctl, group = i, color = burned, fill = burned),
alpha = 0.01) +
geom_density(data = precip.chirps.ww.event %>%
filter(days_before==0) %>%
right_join(slide.regions, by = 'OBJECTID')%>%
inner_join(modis, by='OBJECTID'),
aes(x = pctl, color = 'day_of', fill = 'day_of'), size = 1) +
geom_label(data = cross_df(list(region = region.plot.order,
burned = c('burned', 'unburned'))) %>%
dplyr::mutate(label = c('h', 'i', 'j', 'k', 'l', 'm', 'n',
'o', 'p', 'q', 'r', 's', 't', 'u')) %>%
dplyr::mutate(region = factor(region, region.plot.order, ordered = T)),
aes(x = .95, y = .05, label = label),
hjust = 'right', vjust = 'bottom') +
scale_fill_manual('',
breaks = c('day_of', 'burned', 'unburned'),
values = c('day_of' = NA, burned.colors),
labels = c('day_of' = 'Day of landslide precipitation',
'burned' = 'Bootstrap samples, burned sites',
'unburned' = 'Bootstrap samples, unburned sites'),
guide = guide_legend(override.aes = list(alpha = 0.5))) +
scale_color_manual('',
breaks = c('day_of', 'burned', 'unburned'),
values = c('day_of' = 'black', burned.colors),
labels = c('day_of' = 'Day of landslide precipitation',
'burned' = 'Bootstrap samples, burned sites',
'unburned' = 'Bootstrap samples, unburned sites'),
guide = guide_legend(override.aes = list(alpha = 0.5))) +
scale_x_continuous(breaks = seq(0, 1, 0.25),
labels = c('0', '0.25', '0.5', '0.75', '1')) +
labs(y = 'Kernel density estimate', x = 'Precipitation percentile') +
theme_bw() +
theme(legend.position = 'bottom',
legend.spacing.x = unit(15, "pt"),
axis.text.y = element_blank(),
axis.text.x = element_text(),
axis.ticks.y = element_blank(),
strip.text.x = element_blank(),
panel.spacing.x = unit(15, 'pt')) +
coord_cartesian(expand = F)
plot.bootstrap.sample.kde
plot.bootstrap.sample.kde.legend <- get_legend(plot.bootstrap.sample.kde)
plot.bootstrap.wilcox.legend <- get_legend(plot.bootstrap.wilcox)
cowplot::plot_grid(cowplot::plot_grid(plot.bootstrap.wilcox +
theme(strip.text = element_blank(),
plot.margin = unit(c(5,5,5,5), 'pt'),
legend.position = 'none',
legend.direction = 'horizontal'),
plot.bootstrap.sample.kde +
theme(plot.margin = unit(c(5,5,5,2), 'pt'),
legend.position = 'none'),
rel_widths = c(11, 4)),
cowplot::plot_grid(plot.bootstrap.wilcox.legend,
plot.bootstrap.sample.kde.legend,
rel_widths = c(3, 7)),
nrow = 2, rel_heights = c(10, 0.5))
ggsave('../figures/precip_bootstrap_wilcox+kde.png', width=13, height=11)
The kernel density estimate of the annual distribution of landslides is used to show shifts in seasonality relative to the day of the year.
slide.season <- slide %>%
dplyr::select(OBJECTID, location_accuracy, season, event_yday) %>%
dplyr::right_join(slide.regions, by='OBJECTID') %>%
dplyr::left_join(modis, by='OBJECTID')
plot.seasonality <- bind_rows(slide.season %>%
dplyr::mutate(event_yday = event_yday - 365),
slide.season,
slide.season %>%
dplyr::mutate(event_yday = event_yday + 365)) %>%
dplyr::group_by(region, burned) %>%
tidyr::drop_na(burned) %>%
tidyr::nest() %>%
dplyr::mutate(kde = map(data,
function (df) density(df$event_yday,
from=-100, to=450,
n = 512, bw = 21)),
kde = map(kde, function (dens) tibble(event_yday = dens$x,
density = dens$y))) %>%
tidyr::unnest(kde) %>%
dplyr::filter(event_yday >= 0, event_yday < 365) %>%
ggplot() +
facet_wrap(region~., ncol = 1, strip.position = 'right') +
geom_area(aes(x = event_yday, y = density, group = burned, fill=burned),
trim = T, alpha = 0.6, color = 'black',
position = position_identity()) +
geom_text(data = tibble(text = c('h', 'i', 'j', 'k', 'l', 'm', 'n'),
hjust = c(-3.5, -10, -10, -4, -10, -2, -3.5),
region = region.plot.order,
x = -2.1*365, y = 0.50) %>%
mutate(region = factor(region, region.plot.order, ordered = T)),
aes(x = 0, y = .0035, label = text, hjust = hjust), size = 7) +
coord_polar() +
scale_fill_burned +
scale_x_continuous(breaks = (season.breaks + 365/8) %% 365,
labels = c('W', 'Sp', 'Su', 'F'),
minor_breaks = season.breaks) +
labs(y = '', x = ' ', fill = 'Season') +
theme_bw() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(angle = 0, vjust = 1),
legend.position = 'bottom',
panel.grid.minor.x = element_line(color = 'black'))
plot.seasonality
Smoothed precipitation frequency is used to show shifts in the annual distribution of landslides relative to the annual precipitation seasonality.
n.years <- 2
calc.prop <- function (col) sum(!(col %in% 0)) / n()
calc.yday.frequency <- function (precip, slide, slide.regions, modis) {
precip %>%
dplyr::mutate(precip_yday = yday(date)) %>%
dplyr::select(OBJECTID, precip_yday, pctl) %>%
dplyr::inner_join(slide %>% dplyr::select(OBJECTID, event_yday),
by='OBJECTID') %>%
dplyr::mutate(ydays_before = event_yday - precip_yday) %>%
dplyr::select(OBJECTID, ydays_before, pctl) %>%
dplyr::mutate(ydays_before = ydays_before %% 365) %>%
dplyr::right_join(slide.regions, by = 'OBJECTID') %>%
dplyr::left_join(modis %>% dplyr::select(OBJECTID, burned),
by='OBJECTID') %>%
dplyr::group_by(region, burned, ydays_before) %>%
dplyr::summarize_at(vars(pctl), list(precip.yday.freq = calc.prop))
}
get.data('precip.chirps.yday.frequency',
calc.yday.frequency,
precip.chirps.ww.all,
slide,
slide.regions,
modis)
calc.frequency <- function (precip, precip.yday.freq,
slide, slide.regions, modis) {
precip.chirps.ww.all %>%
dplyr::transmute(OBJECTID, date = lubridate::as_date(date), pctl) %>%
dplyr::inner_join(slide %>% dplyr::select(OBJECTID, event_date),
by='OBJECTID') %>%
dplyr::mutate(days_before = as.integer(event_date - date)) %>%
dplyr::select(OBJECTID, days_before, pctl) %>%
dplyr::filter(days_before < n.years * 365,
days_before > -n.years * 365) %>%x
dplyr::right_join(slide.regions, by = 'OBJECTID') %>%
dplyr::left_join(modis %>% dplyr::select(OBJECTID, burned),
by='OBJECTID') %>%
dplyr::group_by(region, burned, days_before) %>%
dplyr::summarize_at(vars(pctl), list(precip.day.freq = calc.prop)) %>%
dplyr::mutate(ydays_before = days_before %% 365) %>%
dplyr::left_join(precip.yday.freq, by=c('region', 'ydays_before', 'burned')) %>%
dplyr::mutate(precip.freq.norm = precip.day.freq / precip.yday.freq,
precip.freq.anomaly = precip.day.freq - precip.yday.freq,
precip.freq.lt.anomaly = precip.day.freq - mean(precip.day.freq)) %>%
dplyr::ungroup()
}
get.data('precip.chirps.frequency', overwrite = F,
calc.frequency,
precip.chirps.ww.all,
precip.chirps.yday.frequency,
slide,
slide.regions,
modis)
plot.frequency <- function (data, column, y.lable, n.years, groupby=c('burned')) {
data %>%
dplyr::group_by_at(groupby) %>%
dplyr::mutate(smooth = rollmean(!!as.name(column), 90, na.pad=TRUE)) %>%
ggplot() +
geom_line(aes(x = days_before, y = !!as.name(column),
group = burned, color = burned, size = 'daily')) +
geom_line(aes(x = days_before, y = smooth, group = burned, color = burned,
size = 'rolling')) +
geom_hline(data = data %>%
dplyr::group_by_at(groupby) %>%
summarize(mean = mean(!!as.name(column))),
aes(yintercept = mean, color=burned), linetype = 'dashed') +
geom_text(data = tibble(text = c('a', 'b', 'c', 'd', 'e', 'f', 'g'),
region = region.plot.order,
x = -1.95*365, y = 0.50),
aes(x = x, y = y, label = text), size = 7,
hjust = 'right') +
scale_size_manual('',
values = c('daily' = 0.1, 'rolling' = 0.75),
labels = c('daily' = 'Daily frequency',
'rolling' = '90-day rolling average')) +
scale_color_manual('',
values = c('burned' = '#d95f02', unburned = '#7570b3'),
labels = c('burned'='Burned', 'unburned'='Unburned'),
guide = guide_legend(override.aes = list(size = 2))) +
scale_x_reverse('Years before the landslide',
breaks = seq(-365 * n.years, 365 * n.years, 365),
labels = function (breaks) breaks / 365,
expand = expansion()) +
labs(y = y.lable) +
theme_bw() +
theme(legend.position = 'bottom')
}
plot.frequency.region <- precip.chirps.frequency %>%
plot.frequency('precip.freq.lt.anomaly',
'Precipitation frequency anomaly - long-term mean',
n.years=n.years,
groupby=c('burned', 'region')) +
facet_wrap(region~., ncol=1, strip.position = 'right')
plot.frequency.region
cowplot::plot_grid(
cowplot::plot_grid(
plot.frequency.region +
theme(strip.text.y = element_blank(),
legend.position = 'none',
plot.margin = unit(c(5, 1, 5, 5), 'pt')),
plot.seasonality +
theme(legend.position = 'none',
plot.margin = unit(c(5, 5, 12, 0), 'pt')),
rel_widths = c(10, 2.5)
),
get_legend(plot.frequency.region),
nrow = 2, rel_heights = c(10, 0.5)
)
ggsave('../figures/precip_frequency_seasonality.png', width=11, height=12)